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tivity in quasi-two-dimensional Hubbard model. The integral equations for the Green's function 
are self-consistently solved by numerical calculation. Solutions for the order parameter, London 
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are compared with the experiments for the cuprate high-temperature superconductors. Numerical 
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I. INTRODUCTION 

The Hubbard model has been considered as the basic 
model to study the mechanism of high-temperature su- 
perconductivity in the cuprates. 1 By this model, the spin- 
fluctuation-exchange between electrons is considered as 
responsible for the mechanism of high-temperature su- 
perconductivity. A number of calculations, taking into 
account of the spin-fluctuation effects, have been devoted 
to investigating the superconducting properties of the 
two-dimensional Hubbard models. 2-13 

It is proved that the spin-fluctuation theory can suc- 
cessfully describe a number of properties, including the 
temperature dependences of the antiferromagnetic cor- 
relation length 9 and the electric resistivity, 14 of the 
cuprates at high temperatures. However, in most of the 
calculations on the Hubbard model, the superconducting 
pairing is treated by the mean-field like approximation. 
Such an approximation is not appropriate because the 
pairing fluctuation is significant in the low-dimensional 
superconducting systems. 15-20 In fact, the pairing fluc- 
tuation can result in new physical consequence. It is 
believed that the pairing fluctuation is relevant to the 
pseudogap phenomena 18,19 ' 21-24 observed in the normal 
state 25,26 as well as in the superconducting state 27 of the 
cuprates. 

One of the approaches treating the pairing fluctu- 
ation is the ladder-diagram approximation, which has 
been developed on the quasi-two-dimensional (Q2D) 
phenomenological model 18 ' 19 and also on the two- 
dimensional Hubbard model. 20 By the ladder-diagram 
approximation, the long-wavelength fluctuation is taken 
as the predominant contribution. It has been shown that 
the pairing fluctuations can result in considerable reduc- 
tion of the transition temperature T c . According to this 
approach, T c vanishes in the absence of inter-layer cou- 
pling. The reason is that the pairing fluctuation is di- 
vergently strong in the two-dimensional system. This is 
consistent with the Mermin-Wagner-Hohenberg (MWH) 
theorem. 28 



In this work, we intend to study the pairing-fluctuation 
effect on the Q2D Hubbard model. In addition to the 
spin- fluctuation-exchange (S-FLEX), we take into ac- 
count of the contribution from the pairing fluctuation 
in the self-energy of the one-particle Green's function. 
With this spin-pairing-fluctuation-exchange (SP-FLEX) 
approximation, we investigate the superconductivity in 
the Q2D Hubbard model. By self-consistently solving the 
integral equations for the Green's function, we calculate 
the order parameter, London penetration depth, density 
of states (DOS), and transition temperature. Some of the 
results are compared with experiments for the cuprate 
high-temperature superconductors. In the meanwhile, 
we also present some numerical techniques in details in 
the appendices, which is necessary for carrying out the 
numerical solution for the Green's function. 



II. FORMALISM 

The Q2D Hubbard model defined on a layered cubic 
lattice is of the following form 29 

H = - *« c L c Jc« + U ^n iT n u - /i^(n iT + mi) 

(ij)," ' < 

(1) 

where iy denotes the hopping energy of electron between 

the lattice sites i and j, c\ a (ci Q ) represents the electron 
creation (annihilation) operator of spin-a at site i, rii a — 
c\ a c-iu , U is the on-site Coulomb interaction, and n is the 
chemical potential. The (ij) sum runs over the nearest- 
neighbor (NN) sites. In the following, we shall assume 
<y = t for the intra-layer NN hopping and £y = t z for the 
interlayer NN hopping. A quasi-two-dimensional system 
is characterized by the condition t z /t <C 1. Throughout 
this paper, we use unites in which h = ks = 1. 
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A. Normal state 

For simplifying the Feynman diagrams for the Green's 
function, we here firstly present the approximation 
scheme for the normal state. The result for the super- 
conducting state can be immediately obtained by adding 
the anomalous Green's function contributions, and will 
be presented in the next subsection. The normal Green's 
function for the electrons is given by 



G(k, Zn) 



(2) 



where z n = i(2n — l)irT with T the temperature of the 
system is the imaginary fermionic Matsubara frequency, 
£k = — 2t{cosk x + cosfc y ) — 2i z cosfc z — /x, and E(k, z n ) 
stands for the electron self-energy. For brevity, occasion- 
ally, we use the generalized momentum k = (k, z n ) in 
this paper. 

Figure 1 shows the approximation scheme for the self- 
energy. The first two diagrams in Fig. 1(a) are the well- 
known S-FLEX approximation. These two diagrams can 
be combined into a single diagram by redefining an effec- 
tive interaction V e g that is the summation of two inter- 
actions given by Figs. 1(b) and 1(c). The expression for 

V eff is 5 ' 7 ' 8 



V eS (q) 



3 U 2 X (q) , 1 U 2 X (q) 



2 1 + U X (q) 2 1 - U X (q) 



-U\{q) (3) 



with 



^) = §E G ( fc + 9)G(fc). 



(4) 



The generalized momentum q stands for (q, Z m ) with 
Z m = i2mirT the bosonic Matsubara frequency. The first 
and second terms in right-hand-side of Eq. (3) come re- 
spectively from the spin and charge fluctuations. The last 
term eliminates a double counting in the second-order 
diagrams. Owing to the predominant spin fluctuation, 
it is named as spin-fluctuation-exchange approximation. 
The Hartree term has been neglected since it is a con- 
stant, which can be absorbed in the chemical potential. 
In the present calculation, we add the third diagram in 
Fig. 1(a) representing the contribution from the pairing 
fluctuation. Apart from two interaction sides, the shaded 
part essentially represents the processes of the electron 
pair propagating. Figure 1(d) gives the ladder-diagram 
approximation for it with the second order term given by 
Fig. 1(c), which describes the basic propagating process. 
The pairing interaction between two electrons of oppo- 
site spins contains two parts, one due to the transverse 
spin fluctuation (TSF) as given by Fig. 1(c), and another 
one the screened Coulomb potential (SCP) given by Fig. 
1(f). In the right-hand-side of diagrammatic equation of 
Fig. 1 (c), the first diagram represents the propagating 
of a pair without changing their spins in the interme- 
diate state since they interact through SCP during the 



process. In the second diagram, the intermediate spin 
configuration is changed because the two electrons inter- 
act through the TSF. The third diagram describes the 
process the two electrons firstly interact through SCP 
and then through the mediation of TSF, with a minus 
factor stemming from once appearance of TSF. The last 
diagram is similar to the third one but with inverse in- 
teraction sequence. For brevity, we have dropped all the 
momentum on these diagrams. The momentum and spins 
attached to the ladder-diagram are illustrated in Fig. 2. 
For the sake of discussion, we here introduce a notation 
L a f3,/3'a'(k 1 q — k;q — k',k') for the ladder-diagram. In 
following, we will show that at T < T c the value of 
the ladder-diagram diverges at the long-wavelength limit, 
q — > 0. Therefore, the pairing fluctuation represented by 
the ladder-diagram gives significant contribution to the 
self-energy. 
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FIG. 1. Approximation scheme for the self-energy, (a) 
Self-energy for the a-spin electrons. The first term comes 
from the coupling of the a-spin electrons with the density 
fluctuation of opposite /3-spin electrons. The second term is 
due to the coupling between transverse spins through their 
fluctuation. The last term represents the contribution from 
the pairing fluctuation, (b) Interaction between a-spin elec- 
trons due to the density fluctuation of /3-spin electrons, (c) In- 
teraction between transverse spins stemming from their fluc- 
tuation, (d) Ladder-diagram approximation to the pairing 
fluctuation, (e) Second order ladder-diagrams, (f) Screened 
Coulomb interaction between electrons of opposite spins. 

To see how the pairing fluctuation takes effect, 
we consider Fig. 1(d) for the case of a pair of 
electrons with opposite spin and opposite momentum 
that is the case of the ladder-diagram at the long- 
wavelength limit. Because of L Q/ 3 i( 3/ Q /(fc, —fc; —fc', fc') = 
—Lp a! 0' a i(—k,k;—k',k'), we thereby can combine the 
last two terms in Fig. 1 (d) with an effective pairing inter- 
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action Vp denned by Fig. 3(b), and obtain an equation 
like to Fig. 3(c) but with 'w' replaced by '=' under the 
ladder-diagram approximation. To solve this equation, 
we expand the effective pairing interaction in terms of a 
complete set of basis functions (f> n , 



V P (k,k r ) = J2^n(k)Mk')- 



(5) 



(a) 



(b) 
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The function </>„ satisfies the eigen-equation [see Fig. 
3(c)] 

^ Vp(fc, k')G(k')G(-k')(b n (k') = \ n cf> n (k) (6) 



where N is the total number of lattice sites, and A„ is 
the eigenvalue. By so doing, we get 



L a /3,i3'a'(k, —k; —k, k) 



^ Kvn<j>n(k)4> n {k') ^ 

1 — X n 



At the transition temperature T c , the largest eigenvalue 
equals unity, by which the eigen-equation (6) then re- 
duces to the gap equation. At the superconducting state, 
the eigen-equation (6) is modified by adding the term of 
the anomalous Green's functions with the largest eigen- 
value being unity unchanged. Therefore, at T < T c , 
L a f) t f)' a '(k,—k;—k',k') is infinitive. It implies that the 
long- wavelength pairing fluctuation gives significant con- 
tribution to the self-energy. On the other hand, as T — > 
T c from the normal state, with the largest eigenvalue of 
Eq. (6) approaching to unity, L a p^> a >{k, — fc; —k', k') di- 
verges at this limit. Therefore, for the normal state at 
T close to T c , the long- wavelength pairing fluctuation is 
important as well. 



q-k'fi 




k,a k'.a 
FIG. 2. Ladder-diagram representing the propagating of 
a pair of total momentum q. k' and q — k' are the initial 
momentum of the a and /3-spin electrons, respectively. After 
propagating, their momentum change to k and q — k. 

From the right-hand side of Eq. (7), we see that ex- 
cept for the term corresponding to the largest eigenvalue, 
all other terms are finite. We therefore keep only the 
most diverging term to simplify the solution of the equa- 
tion given by Fig. 1(d). That is, we can consider only 
the pairing of largest eigenvalue. 24 For the present case, 
the largest one is the d-wave pairing. Hereafter, we de- 
note the largest eigenvalue and the corresponding eigen- 
function simply as Xd and 4>{k), respectively, and the 
coupling constant simply by v. 



(C) 



FIG. 3. Approximation to the ladder-diagram, (a) Approx- 
imate second order ladder-diagram obtained from Fig. 1(e) 
by neglecting the dependence of the interactions on the total 
momentum of the pair, (b) Effective pairing interaction, (c) 
Renormalized diagrammatic equation for the ladder-diagram. 

On observing the above-mentioned fact, we make the 
approximation in Fig. 1(d) using the pairing inter- 
actions of zero total momentum because near q = 
the pairing fluctuation is most significant. We then 
obtain equations as given by Fig. 3 for determining 
L a p,f3' a '{k,q — k;q — k',k'). Further more, by consid- 
ering only the d-w&ve pairing that has the largest eigen- 
value, the ladder-diagram given by Fig. 3(c) reduces to 
the same one as we previously encountered for the phe- 
nomenological model. 19 Applying the previous result to 
the present case (see Appendix A), the ladder-diagram 
summation is obtained as 



Lal3,p a {k, q-k;q-k', k') = P(q)(j>(k)4>(k'), 



with 



v 2 n( q ) 



v 2 n(q) 



1 + vU(q) 



(8) 

(9) 
(10) 



In Eq. (9), the last term eliminates all the second-order 
diagrams since the contribution to the self-energy from 
the first two diagrams in the right-hand-side of the di- 
agrammatic equation of Fig. 1(e) has been taken into 
account in S-FLEX approximation, while the other two 
diagrams are negligible with comparing to the infinitive 
ladder-diagram summation. For the self-energy, the final 
expression is 

s ( fc ) ^ -jj E G ( fc - ^ ff (<z) + E G (« - k ) p M 



(ii) 



As we have noted above, the last term in Eq. (11) is cor- 
responding to the previous approximation 19 for the phe- 
nomenological model. 18:19 In that model, however, the 
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interaction is simply a constant ci-wave pairing potential, 
and the hopping energy is proportional to the hole con- 
centration for taking into account of the constraint ex- 
cluding the double occupation. The prohibition of dou- 
ble occupation stems from the t — J model that is the 
large U limit of the Hubbard model. It is a consequence 
of strong short-range antiferromagnetic coupling. Un- 
der the SP-FLEX approximation scheme for the Hubbard 
model, however, the antiferromagnetic coupling is taken 
into account by the first term in Eq. (11). Moreover, for 
not too large U, the d-w&ve pairing potential given above 
varies with the temperature, hole concentration and U. 

By passing this subsection, we give an expression for 
the pairing potential. As seen from Figs. 1(c), 1(f) and 
3(b), Vp(k, k 1 ) equals V 1(c) {k + k')-V m {k-k'). Since 
the pairing function (f>{k) can be taken as even function of 
k, the pairing potential in Eq. (6) can be thereby written 
as Vp(k — k'). We then have 



1 U 2 X (q) 



2 1 + U X (q) 2 1 - Ux(q) 



u. 



(12) 



Therefore, the left-hand-side of Eq. (6) is a convolution 
of Vp and the rest part. 

In addition, the chemical potential \x should be deter- 
mined to yield the hole concentration 



5 = -^E[ G w + G (- fc )]- 



(13) 



All the above equations form the closed system that self- 
consistently determines the Green's function. 



B. Superconducting state 

For the superconducting state, the above results should 
be extended to including the contributions from the 
anomalous Green's function. In the Nambu represen- 
tation, the Green's function is given by 



G(k) = [Z n - - t(k)Y 



(14) 



where z n is understood as z n cro, and a is the Pauli ma- 
trix. Occasionally, we will use the Pauli components of G 
defined by G = Go + Gi<7i + G 3 a 3 . Correspondingly, the 
self-energy is expressed as £ = So + Si<7i + £303 as well. 
The diagonal element En(fc) = S (fc) + S 3 (fc) is given 
by the same diagram Fig. 1(a) except where the effec- 
tive interaction and the ladder-diagram should include 
the contribution from the anomalous Green's function. 
The element £22^) is obtained by E 2 2(fc) = -En(-fc). 
The off-diagonal part Si is given by the gap equation 



(15) 



The expressions for V e s and Vp can be obtained as 



V eB (q) = 



3 U 2 X -(q) , 1 U 2 X+ (q) 



21 + U X -(q) 2 1 -U X +(q) 



U 2 X (q) (16) 



Vp(q) = 



U 2 X -(q) 1 U 2 X+ (q) 



21 + u X -(q) 21 -U X +(q) 



+ U 2 X l(q)-U 



with x±(q) = x(q) ± Xi(q), and 



(17) 



(18) 



The expression for xil) is the same as Eq. (4) where the 
Green's function G(k) is understood as Gn(fc). A simple 
derivation of these results is presented in Appendix A. 

Following the similar analysis as in the previous subsec- 
tion, one can take a corresponding approximation for the 
pairing fluctuation. In the superconducting case, how- 
ever, besides the diagonal pair propagating (pairing of 
particles or holes), attention must be paid to the off di- 
agonal pair propagating as well. The latter is the process 
that the initial state is a pair of particles (holes), while 
the final state is a pair of holes (particles). Therefore, the 
pair propagators satisfy a matrix equation. The ladder- 
diagram approximation with ci-wave channel interaction 
is given in Appendix A. The function P(q) appeared in 
En represents the pair propagating. It can be divided 
into two parts, P{q) = Po(q) + Ps(q)- Their expressions 
are given by 



Po(q) 

Ps(q) 



v[D(q) - 1 
v 2 U 3 (q)[l- 



~vU a (q)}/D(q)-v 2 U Q (q) 
D(q)]/D(q) 



D(q) = [1 + vU+(q)][l + vU-(q)] - v^q) 
with n±(g) = U (q) ± Ili(g), and 



(19) 
(20) 
(21) 



T 



n °(«) = n E ^( fc )[Go(fc)G (fc -q)- G 3 (k)G 3 (k - q)} (22) 



N 
T 



(23) 



U ^ = n E <t> 2 {k)[G 3 (k)G (k -q)- G (k)G 3 (k - q)]. (24) 



The eigen-equation for determining the function <f>(k) 
now is extended to 



T 

N 



V P (k k')[G 2 (k') + G\{k') - Gl{k')]cj){k') = m- 

k' 

(25) 



with the largest eigenvalue being unity. This equation is 
equivalent to Eq. (15) since 4>(k) differs from Si(fc) by a 
normalization constant. It also leads to 1 + t>n_(0) = 0, 
and thereby P(q) diverges at q = 0, which means the 
existence of the Goldstone mode. 
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In terms of the above functions, the diagonal parts of 
the self-energy can be expressed as 

1 

_T^kl^ [Go{k q)Po(q) ^ qWg)l 
q 

S 3(fc)--^E G 3( fc -^ff(9) 



l + vn ± (q) 



l + vn±(q)+cfq 2 z 



with 



(26) 



n±(g) = n±( g )|, 

_i_ v d 



W n±(<Z) '- 



(31) 

(32) 
(33) 



N 
<i 

T(f> 2 (k) 
N 



and use fl±(q) for Ii±(q) in the nominator in Eq. (30). 
Note that cf is defined as the derivation in the right-hand 
side of Eq. (33) at q = since where the q z dependence 



y^[Go(fc — q)P3{q) ~ G 3 (k — q)Po(q)] (2i^)i m portant. To evaluate the constants cj, we need to 



By using the Pauli component of the Green's function, 
the expression for Eq. (13) can be simplified as 



(28) 



So far, we have all the equations for the superconducting 
case. 



C. Q2D approximation 

The Green's function G(k) and the susceptibilities 
x(<?)'s and n(q)'s are defined in three-dimensional space. 
Actually, in case of t z /t <C 1, they very weakly depend 
on the ^-component variables. However, the dependence 
on q z of function P(q) is delicate. Consider the denomi- 
nator function D(q) at small q and Z m = 0. Since 11(g) 's 
are even functions of q, we have 



D(q,0)nc(q 2 x +ql) + c z q 2 



(29) 



where c and c z are constants. The g^-term in Eq. (29) 
comes from the interlayer electron hopping. The ratio 
c z /c is much less than unity. If the (j z -dependence in 
D(q) is ignored, the second summations in Eqs. (26) 
and (27) will be divergent, which implies there will be no 
superconductivity in the system at finite temperature. 
This conclusion is consistent with the MWH theorem. 28 
We therefore need to keep the ^-dependence in the de- 
nominators of P(q)'s at least to the order q 2 . 

For illustrating our approximation scheme, we firstly 
consider the case of Z m = 0. Since n 3 (q, 0) = 0, we have 
P 3 (q,0) = and 



Pn = 



IL 



+ 



IE 



2 L l + vIL- 1 + vU+ 



v 2 U 



(30) 



where the arguments (q, 0) have been dropped for 
brevity. As has been mentioned in last subsection, the 
first denominator l+uTI_ vanishes at q = 0. Even though 
the second denominator 1 + vll + is finite at T < T c , it is 
small. Especially, it vanishes too at T — T c . Therefore, 
we expand both of the denominators to the order q 2 , 



take the derivative of the Green's functions G{q — k)'s 
with respect to q z as indicated by Eqs. (22) and (23). By 
neglecting the q z -dependence of the self-energy, G(q—k)'s 
thereby depend on q z only via It is expectable 

that such an approximation does not change the physical 
result so much. To the second order of t z /t, we obtain 19 



tjvT 
N 



k 



with 



<?{k){[^-G 3 {k)] 2 T [^-Gmf ~ bJ-G (fc)] 2 }, 

OC,k O^k d£k 

(34) 



— Go(fc) « 2G (fc)G 3 (fc) 

O^k 



d_ 

Wk 

d_ 

d 
Wk 



^-Gi(fe) « 2Gi(fc)G 3 (fc) 



G 3 (k)^G 2 (k)-G 2 (k) + G 2 (k). 



(35) 



With such an approximated Po(q, 0), the integral over q z 
in Eqs. (26) and (27) at Z m — can be taken imme- 
diately by neglecting the ^-dependence in the Green's 
function. Therefore, instead of Po(q, 0) in Eqs. (26) and 
(27), we insert in a function defined by 



P o cff (q,0) = - r^P (q,0) 

71" JO 



= y [n_/_ + n + / + ]- w 2 n , (36) 



with 



f± 



7 ± (5) 



7TC Z 



atan[7T7 ± (g)] 



ll/2 



i + «n±( g ) J 



(37) 
(38) 



where again, in the last line of Eq. (36), we have dropped 
the arguments (q, 0) for brevity. 

We now consider the situation of Z m ^ 0. D{q) is 
finite in this case. However, in consistent with the ap- 
proximation for P (q, 0) , we still keep a small q^-term 
in the denominator D{q). Though this g z -dependence is 
negligible at high temperature, it is reasonable in case 
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of low temperature. According to the expansion by Eq. 
(31), we expand D(q) as 

D(q)=D(q)+e(q)q 2 (39) 



with 

e(q) = [1 + vTL + (q)]c- + [1 + vIL_ (g)]c+. (40) 

To the order q 2 zl this expansion reduces to the result for 
the case of Z m = 0. Correspondingly, we can define the 
functions Pg ff (<7) and P| ff (q) by taking the integral over 
q z . This procedure is equivalent to replacing 1/D(q) in 
Eqs. (22)- (24) with a function f(q) defined by 



/(</) 



7(9) 



atan[7T7((j')] 
e(<?) 



?re(g) 



L ^) J 



(41) 
(42) 



With the functions P(<z)'s in Eqs. (26) and (27) replaced 
with P eff (g)'s, the problem of numerically solving the 
integral equations is then reduced to a two-dimensional 
one. All discussed above are for the case of supercon- 
ducting state. The results for the normal state can be 
obtained by setting G\(k) = 0. 



mesh in momentum space for the inverse transform. The 
values of V e s{<\, 0) for this mesh are obtained by local 
quadratic polynomial interpolation of the smooth func- 
tions x(<z)' s given on a 128x128 mesh. On the other 
hand, the function P§ a (q, 0) has divergently sharp peak 
at q = and T < T c . In Appendix C, we deal with the 
inverse transform of this function. 

The difficulty in solving the eigenvalue problem given 
by Eq. (6) is that the memory requirement for the coef- 
ficient matrix is huge. It is impossible to solve this equa- 
tion in momentum space. In Appendix D, we rewrite 
the eigen-equation in real space. At high temperature 
not too close to T c , Eq. (6) can be solved in real space 
with a small number of lattice sites in a reduced region. 
This reduces greatly the amount of numerical calculation 
work. 

The integral equations determining the Green's func- 
tions are numerically solved by iteration method. Once a 
solution at temperature T is obtained, it is then used as 
an initial input for the next calculation at temperature 
T + ST. A more efficient way is to use an extrapolation 
from the known solutions at temperatures T\ and T 2 as 
the initial input for the next solution at T2 + ST. 

In the present calculation, we set U/t — 5 and t z /t — 
0.01. All the results presented in the figures are for these 
parameters. 



III. NUMERICAL RESULTS 

Since the functions G(k), V c g(q) and P cff (g)'s are de- 
fined in multi-dimensional space, they require huge mem- 
ory storage in the numerical computation process. Espe- 
cially, the function PQ S (q) is singular at q = 0. There- 
fore, to carry out the numerical solution, we need to de- 
velop numerical scheme to reduce the amount of compu- 
tation without losing the accuracy. In Appendix B, we 
present our scheme for the Matsubara frequency summa- 
tion. The summation is taken over 57 points, a subset 
of the frequencies, in a sufficiently large range. The cut- 
off frequencies are z c = (2N C — l)nT for the fermions 
and Z c = 2(N C — l)nT for the bosons, respectively, with 
N c = 1017. For the typical temperature T/t ~ 0.01 un- 
der consideration, this means 2N c nT/t <~ 64. For calcu- 
lating the function x{l)i beyond this range, the summa- 
tion over the terms of n > N c is analytically carried out 
by using the asymptotic formula of the Green's function 

G(k, z n ) -» l/z n . (43) 

The error of the summation over the terms of n > N c is 
of the order 0(z~ 3 ). 

The convolutions in the momentum space are carried 
out with fast Fourier transforms (FFTs) on a 128x128 
lattice. For inverse transform of T4ff(q, 0) and P§ R (q, 0), 
i.e., from momentum space to real space, we have to pay 
special care. At low temperature, V e s(q, 0) has strong 
peaks near q = (it, ir). 7 We therefore use a 256x 256 



A. Eigenvalue Ad 

In Fig. 4, we show the result for the eigenvalue A<j as 
function of T at S = 0.125. The S-FLEX result is also 
presented for comparison. Due to the pairing fluctuation, 
the eigenvalue by the present SP-FLEX approximation 
is considerably reduced from that of the S-FLEX, giving 
rise to a lower transition temperature. Moreover, there 
is a distinguishable difference between their behaviors at 
temperatures close to T c . By the S-FLEX approximation, 
we have d\ d (T)/dT ^ at T = T c . It means that at 
T < T c , by keeping no superconducting pairing, the S- 
FLEX approximation allows a solution of A^ > 1. In 
contrast to this feature of the S-FLEX approximation, 
the solution for A<j by the SP-FLEX approximation can 
never get across the line A<2 = 1. At T — > T c , the pairing 
fluctuation effect is more pronouncedly with A<j — > 1, 
which conversely suppresses Ad. As a result, the curve 
Ad(T) is smoothly connected to the straight line A^ = 1. 
That is 

^A d (T)| Tc =0. (44) 



The inset of Fig. 4 shows that y/l — Xd varies nearly 
linearly as T — > T c , which means 1 — A^ oc (T — T c ) 2 . 

A problem then comes in the determination of T c by 
A(T C ) = 1 from the side of normal state. Because of 
dT/dXd = 00 at T = T c , a small numerical error in A^ 
may results in considerable error in T c . Therefore, T c 
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cannot be accurately determined by the function A<j(T). 
In our numerical calculations, we solved Eq. (6) using 
two sets of numbers, 25 and 49 respectively, of the lat- 
tice sites in the reduced region (see Appendix D). The 
solid circles in Fig. 4 represent the numerical results of 
49 lattice sites, while the circles are the ones of 25 lattice 
sites. Close to T c , the difference between the two results 
are visible. Even with 49 lattice sites, the transition tem- 
perature so determined is not reliable. 
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FIG. 4. Eigenvalue Ad as function of T at U/t = 5 and 
S = 0.125. SP-FLEX is the present approximation. The result 
of S-FLEX is also depicted for comparison. The solid circles 
represent the numerical solution to Eq. (6) solved with 49 
lattice sites in a reduced region (see Appendix D). The circles 
are results of 25 sites. The squares are the transition points 
obtained by Eq. (25). The inset shows (1 — Ad) 1 / 2 of the 
SP-FLEX approximation as function of T. 

The problem can be resolved from the superconduct- 
ing side. Instead of Eq. (6), we solve Eq. (25). Since 
the eigenvalue A^ = 1 is known at T < T c , Eq. (25) can 
be solved via iteration on the whole lattice. The squares 
in Fig. 4 denote T c 's obtained from the superconducting 
side for SP-FLEX and S-FLEX, respectively. The tran- 
sition temperatures so determined should be reliable. 



B. Order parameter 

At low temperatures, we have obtained self-consistent 
solutions in which Ei(fe) is finite. The symmetry of pair- 
ing is d-wave, Si(k x , k y , z n ) = -Ei(k y , k x , z n ). Here, 
we define the order parameter, 



with 



A = Ei(X,*i)/Z(X,*i) 



Z{k,z n ) = 1 - E (k, z n )/z r , 



(45) 
(46) 



and X = (jr, 0). The quantity A is a measure of the su- 
perconducting gap at the Fermi surface near the point 
X. 7 ' 8 Fig. 5 shows the order parameter as function of T 
at various doping concentrations. The symbols represent 
the numerical data, while the lines are the extrapola- 
tions. As seen from Fig. 5, close to T c , A decreases dra- 
matically. (In the numerical calculation, because of this 
rapid decreasing, to get an iteration converged at tem- 
perature T + ST with an initial input extrapolated from 
some other solutions at and close to temperature T, the 
change ST must be very small. At low doping concen- 
trations, close to T c , a change of ST jT <~ 10~ 4 at most 
was allowable in the present calculation.) At T = T c , we 
have dA(T)/dT = oo. Because of this divergence, T c can 
be accuratclv determined bv A = 0. 
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FIG. 5. Order parameter A as function of temperature T 
at various hole concentrations S. The symbols represent nu- 
merical data. The lines are extrapolations. 
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FIG. 6. Ratio 2A / k B T c at various hole concentrations <5. 
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On the other hand, at low temperature, A(T) should 
be flat. We then can obtain the value Ao = A(0) by 
extrapolation. In Fig. 6, we show the ratio 2Ao//cbT c 
at various doping concentrations. At the undcrdoping 
regime, the ratio is about order 10. It decreases with 
the doping concentration. This ratio is considered as a 
characterization of the coupling strength of the system. 
The low doping regime corresponds to strong coupling. 
With increasing the doping concentration, the coupling 
gets becoming weak. 



C. Phase diagram 

As mentioned in the previous subsection, the transi- 
tion temperature T c is determined by A(T C ) = 0, which 
gives accrate solution for T c . The present SP-FLEX re- 
sult at U/t = 5 (solid circles) for the boundary of su- 
perconducting phase in the T c — S phase diagram is de- 
picted in Fig. 7. All the solid circles are obtained by 
numerical solutions. The S-FLEX result (circles) is also 
shown for comparison. Clearly, due to the pairing fluc- 
tuation, T c is considerably reduced from that of the S- 
FLEX approximation. The reduction is more significant 
at lower hole doping where T c decreases with decreasing 
5. On the other hand, at large 5, the pairing fluctuation 
is less pronouncedly. This is consistent with the previ- 
ous conclusion. 19 The results of the previous calculation 
on the phenomcnological model (with coupling constant 
J/t = 0.2) 19 and the experiments 30 for the cuprate high- 
temperature superconductors are exhibited in the inset 
of Fig. 7 for comparison. The behavior of the previous 
result at small hole doping clearly differs from the ex- 
periment and the present calculation. It may stem from 
the crude treatment of the short-range antifcrromagnctic 
coupling by the phenomenological model. In contrast to 
the previous result, the present calculation gives a reason- 
able description of the experiment at small hole doping. 

The dashed line in Fig. 7 is an extrapolation of the 
numerical result. Unlike the previous case for the phe- 
nomenological model, 19 the extrapolation gives a nonzero 
minimum hole concentration, S m very close to 0.05, of the 
phase boundary. Again, this seems reasonably reflecting 
the feature of the experimental result. The largest transi- 
tion temperature T c max obtained by SP-FLEX is about 
0.0137i at 5 = 0.175. Using t w 0.6 eV, 18 < 19 we have 

T RjQSX 
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FIG. 7. Transition temperature T c as function of hole con- 
centration 8. The solid-circles and circles are obtained at 
U/t = 5 by the SP-FLEX and S-FLEX calculations, respec- 
tively. The dashed line is an extrapolation of the numerical 
result. In the inset, the symbols denote the experimental 
results for the cuprate high-temperature superconductors, 30 
while the solid line is the previous result 19 for the phenomeno- 
logical model. 

To see why the pairing fluctuation results in reduc- 
tion of T c , we analyze the self-energy given by Eq. (11) 
[which is the same as that given by Eqs. (26) and (27) 
at T = T c ] with P(q) replaced by P ctf (q). At Z m ~ 0, 
since V e s(q) and P eB (q) have strong negative peaks re- 
spectively at q w Q = (71", 7r) and q = 0, we here make a 
crude approximation (for the sake of illustration) for the 
self-energy, 

S(fe)w -G(k-Q,*„)[& £ U cff (<z)] 



iG(k rZ „)[^Ef cff b)]- 



(47) 



q~0 



The summations in the 
two negative quantities. 
G(k-Q,z„) « -G(k,- 
these facts, we get S(fc) 



square brackets give rise to 
At small 5, fj, w 0, we have 
z n ). Taking into account of 

« -G(k,-z„)r| with r| = 



a\ + ci2(/) 2 (k), where oi and ai are two positive constants. 
Substituting the result into Eq. (2), one obtains 



G(k, Zn) 



211 



(1-W1 + 



_2 >■ 



(48) 



Applying these results to Eq. (6), we see that the factor 
G(k)G(—k) is reduced, and so is T c . 
Physically the quantity 



N 



g~0 
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is a measure of the density of pairs at their excited states. 
Since the function P (q) is given in terms of fi-(q), we 
analyze fi(q) — H(q)\ q ^ = o especially at q = 0. From Eq. 
(10), we see that the quantity 11(0) comes predominately 
from the summation over the points close to the Fermi 
surface in momentum space. At smaller 5, the Fermi 
surface of the Hubbard model is larger, so is the number 
of the pairs at their fluctuating states. Therefore, the 
reduction on T c is larger at smaller 5. 

The form of the Green's function given by Eq. (48) 
implies that there exists a pseudogap in the energy spec- 
trum of the electrons. 19 ' 31 Consider the spectral function 

AQl,E)= -±ImG(k,£ + 



which is nonzero only for E 2 - 4r^ < £ 2 < E 2 . The 
noninteraction delta-function peak becomes a square root 
singularity. Because of the constraint, the area of the k- 
space of A(k, E) ^ decreases at E — > 0, resulting in 
a suppression of the density of states. Especially, the 
density of states vanishes at E = (for /i = 0) since 
the area becomes to zero and the singularity disappears 
there. 

On observing the function T 2 , we note that the pseu- 
dogap stems from the spin and pairing fluctuations. Even 
at T = 0, the pseudogap remains in the superconduct- 
ing state because the spin fluctuation (owing to which 
the superconductivity takes place) and the pairing fluc- 
tuation (coming from the Goldstone mode 19 ) exist. This 
may explain the recent experiment. 27 



D. Density of states 

The density of states is defined by 

k 

The Green's function needs to be analytically continued 
from the imaginary Matsubara frequency to the real fre- 
quency. In terms of an effective self-energy £(fc) defined 
by 



t(k) =S (fc) + S 3 (fc) 



£ 2 (fc) 



z n - 6c - So(fc) + X 3 (fc) ' 
the Green's function Gn(fc) is written as 

Gn(fc) = 



(50) 



(51) 



Using Pade approximation, 32 we have obtained the ana- 
lytical continuation for the effective self energy E(fc). 
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FIG. 8. Density of states p{E) at 8 = 0.125 at various 
temperatures. For clarity, the y-axes for the results of 0.7T C 
and 0.4T C have been offset upwards 0.05 and 0.1, respectively. 
For 1.01T C and 2T C , the offsets are downwards 0.05 and 0.1, 
respectively. 

The results for the density of state at S = 0.125 at 
various temperatures are depicted in Fig. 8. At low tem- 
perature, the width of the gap is nearly constant. Be- 
low T c , a peak-dip-hump (PDH) structure is clearly seen 
above the Fermi energy. The positions of the peak and 
hump are about A and 3 A, respectively. Such a phe- 
nomenon has been observed in the cuprates experiments 
and has been explained by model calculations. 33 ' 34 The 
PDH stems from a coupling between the electrons and a 
collective mode of energy about 2A. Below T c , the super- 
conducting gap opens with the maximum of 2A appeared 
in the region near the points (±7r, 0) and (0,±7r) in the 
Brillouin zone. The collective spin fluctuation mode of 
energy 2A and momentum Q can then exist in the sys- 
tem. An electron of energy 3A can transit to a state 
of energy A (at which the DOS has a peak) by excit- 
ing the collective mode and losing energy 2A. Effec- 
tively, the lifetime of the electrons of energy 3 A is short, 
and thereby the dip appears in the DOS. In the FLEX 
scheme, such a collective mode is described by the effec- 
tive interaction V e g that is self-consistently determined 
by the present calculation. 

On the other hand, above T c , there still remains a pseu- 
dogap in the DOS. With increasing the temperature, the 
minimum moves to high energy. As stated in last sub- 
section, the pseudogap comes from both of the spin and 
pairing fluctuations. At higher temperature, the pairing 
fluctuation is less important. The electron states k and 
k + Q couple with each other mainly through the spin 
fluctuations. Especially, the degeneracy of those states 
at the magnetic zone surface is lifted and the weights 
shift away, resulting in a reduction in the DOS at the 
corresponding energy. 



9 



E. London penetration depth 



In this subsection, we study the magnetic penetration 
depth. The London penetration depth Xl in x-direction 
is given via 

. 2 Anne 2 47r , 

A L 2 = ^ + ^C%)|^o (52) 

with C(g) the y-component current-current correlation 
function defined by 

C(q,T-r') = - <T T [J„(q,r)J !/ (-q,r')] >, (53) 

where < • • • > means a statistical average, T r is the 
imaginary time-r ordering operator, and J y is the y- 
component current operator, 

J(q) = -e^V k ek4_ q/2a c k+q/2a . 

ka 

In Eq. (52), n and m* are the number density and the 
effective mass of electrons, respectively. By using the 
Ward identity for the current vertex, the quantity Co = 
C(q)\ q= o can be written as 

C 'o = ^Et G n( fc ) + Gi2(*)]V k & • v fc . (54) 
k 

with Vfc = Vk[£k + Sn(fc)]. By noting the following 
equations, 

G? 1 (fc)V k [Ck + E(fc)] = VkGn(fc), 

^E V ^-V k G n (fc) = — ^, 
TV ^ to* 

we get the expression for A^ 2 , 

A Z 2 = Et G ? 2 ( fc K - G?i(fc)u fc ] • V k &, (55) 

fe 

with 

= Sf(fc) 
fc k z n -e k -S (fc) + S 3 (A : )- 

By Eq. (55), A^ 2 vanishes identically at and above T c . 




FIG. 9. Quantity \ 2 L (0)/ \ 2 L (T) as function of temperature 
T at various doping concentrations S. The solid circles de- 
note the experimental results for a-axis penetration depth of 
YBa 2 Cu 3 6 .95- 

Shown in Fig. 9 are the results for the quantity A£ 2 (T) 
as function of temperature T at various doping concen- 
trations. The experimental result for the a-axis pen- 
etration depth of YBa 2 Cu306.95 is also presented for 
comparison. 35 It is well known that A£ 2 (T) under d-wave 
pairing symmetry varies linearly with T at low T. Using 
this property, we infer the value A^ 2 (0) from extrapola- 
tion of the known results at finite temperatures. We then 
get the zero-temperature superfluid density n s . Figure 10 
shows the relationship between T c and n s at a number of 
doping concentrations. This result resembles the experi- 
mental observation by Uemura et at. 36 They have found 
a universal linear relation between T c and n s /m* at un- 
derdoping regime. This behavior cannot be explained by 
the BCS theory, nor by the S-FLEX scheme in which T c 
and n s do not decrease with 5 decreasing. 

In the mean-field theory, the electron pairs are con- 
sidered as all in the Bose- Einstein condensate below T c . 
Since the pairing interaction (stemming from the spin- 
fluctuation exchange) is stronger at smaller doping con- 
centration, more electron pairs are produced in the con- 
densate. This leads to higher T c and larger superfluid 
density. In contrast to the mean-field theory, in the SP- 
FLEX scheme, the electron pairs are allowed to occupy 
their excited states. At low temperature, those collective 
modes are the most available excited states. Even at the 
ground state, there remains the zero-point motion for the 
collective modes. So, only a part of the pairs stay in the 
condensate. As stated earlier, the pairing fluctuation is 
stronger at smaller S, resulting in lower T c and lower n s . 
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FIG. 10. T c vs. n 8 at various doping concentrations 8. 
IV. SUMMARY 



Considering the convenience for readers, we here 
present simple derivations of the effective interaction V e ff, 
the pairing potential Vp and the ladder-diagram approx- 
imation for the pair propagators for the superconducting 
state. 

Firstly, we consider the extension of Fig. 1(b). The 
interaction 14 (<?) between a-spin electrons can be written 
as 



(Al) 



where X^^{o) is the density-density response function be- 
tween the opposite /3-spin electrons. Generally, the func- 
tion x aa in the imaginary time-r space is defined as 

X aa ' (q, T - t') = - < T T n a (q, r)n a , (-q, r') > /N 

(A2) 

where n Q (q, r) is the density operator of a-spin electrons. 
In the Mutsubara-frequency space, the Dyson equation 
for x aa (q) reads 

f°' (q) = X aa ' (q) + E X al (q)Ux- W (q) (AS) 



In summary, we have investigated the d-wave super- 
conductivity in the quasi-two-dimensional repulsive Hub- 
bard model. Both of the spin and pairing fluctuations 
are taken into account in the self-energy. We have self- 
consistently solved the integral equations for the Green's 
function. The present calculation reflects a number of 
features of the experiment results for the cuprate high- 
temperature superconductors. The calculated boundary 
of superconducting phase shows a parabolic-like shape, 
reasonably describing the experiments. The peak-dip- 
hump structure in the density of states is naturally re- 
produced by the present calculation. In addition, the 
present calculation gives reasonable explanations for the 
temperature dependence of the penetration depth and 
the relationship between T c and the supcrfluid density. 

The pairing fluctuation implies that an amount of pairs 
occupy their excited states. This fluctuation effect leads 
to the reduction of the condensation of the pairs and 
thereby the transition temperature. On the other hand, 
the spin fluctuation produces the attraction between the 
electrons, resulting in the electron pairing and conden- 
sation. In the meanwhile, the spin and pairing fluctua- 
tions commonly bring about a pseudogap to the electrons 
whether the temperature is below or above T c . 
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APPENDIX A 



where the 7-summation runs over the up and down spins, 
and x's are the irreducible response functions. Under the 
ring-approximation, from Eq. (A2), x QQ ( < ?) is obtained as 
the same as given by Eq. (4) where the Green's function 
G(k) is understood as Gn(fc). x /3/3 ( < z) is equal to x aa (q) 
because of the equality of the up and down spin electrons. 
For a 7^ (3, we have X (o) = X (q) = Xi(<z) as given by 
Eq. (18). Solving Eq. (A3) and substituting the result 
into Eq. (Al), we have 



1 



u x+ (q) i + u X -(q) 



(AA) 



with x±(q) = x(q) ± xi(<?)- 

Secondly, for the interaction between transverse spins 
corresponding to Fig. 1(c), we write 



v c (q) = u 2 x + -(q) 



(A5) 



with x + (q) the response function between transverse 
spins defined by 

X+-(q,T-r') = - < T T S+(q,T)<?(q,T') > /N (A6) 
and S* + (q, r) = J2 c kT( T ) Ck +ql( T )- F rom Eq. (A6), the 

k 

irreducible function can be obtained as x^ (q) = X-{q)- 
Since X + ~(<z) satisfies the Dyson equation x H (q) = 

X-(q) - X-(q)Ux + ~(q), we have 



U 2 X-(q) 
l + U X -(q)' 



(A7) 



With the above results, for the effective interaction 
V eS (q) = V h (q) + V c (q) - U 2 x{q), we get the expression 
as given by Eq. (16). 
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Analogously, the interaction corresponding to Fig. 1(f) 
can be written as 



V l (q) = U + U 2 x aP (q)- 
From Eq. (A3), we have 

ir x+(q) x-(q) 



x af3 (q) = 



2 l i-u x +{q) i + u x -{q) 



(AS) 



(A9) 



The result for the pairing potential Vp(q) = V c (q) — 
Vf(q) + U 2 Xi(q) is then obtained as given by Eq. (17). 
Again, the term U 2 xi(q) eliminates a second-order dou- 
ble counting. 

Finally, the function P(q) appeared in the element 
En(fc) of the self-energy can be expressed as 



P(q) = v 2 [fl n (q)-U n (q)} 



(AW) 



where Ilii(g) is the pair-pair response function, and 
Ilii(g) is the irreducible part (or the pair susceptibility, 
which eliminates the second-order double counting) . The 
function n M1/ (q, t — t') is defined as 

n„„(q,r-T / ) = - < r T ^(q,r)pt( q , r ') > /jV, (All) 

with n,v = 1, 2, pi(q) = £ </>(fc)c q _ kJ c kT and p 2 (q) = 

k 

^2 ^(^) c k-qT c -ki- Strictly speaking, p\ and p 2 are not 

k 

the Schrodinger operators since (p(k) depends on the Mat- 
subara frequency. Here, we formally regard them as only 
depending on the momentum. At the end, we extend 
the result to include the frequency dependence. (Alter- 
natively, one can draw the Feyman diagram from the be- 
ginning. The final result is the same.) From Eq. (All), 
the expressions for the irreducible susceptibilities can be 
obtained as Eqs. (22)-(24). The Dyson equation for 
in matrix form is 



IL(q) = U(q) - vtl(q)fl(q). 



(A12) 



The diagonal parts of the Pauli components of the matrix 

P(q) = v 2 [Il(q) — U(q)] can be expressed as Eqs. (19) and 
(20). 

APPENDIX B 

In this appendix, we intend to develop an algorithm for 
the approximate summation of a series. It is analogous 
to Simpson's integral method. Firstly, we consider the 
following summation, 



5(n ,n 2 ) = 



(Bl) 



where n 2 = n + 2h with h an integer. Suppose f(x) is 
a smooth function over the range no < x < n 2 . We then 
expand f(n) as 



where c\ and c 2 are constants. With the given values 
fj = f(n +jh), for j = 0, 1 and 2, the constants can be 
expressed as 

ci - (-3/ + 4/i-/ 2 )/2/i 
c 2 = (/o - 2/i + f 2 )/2h 2 . 
Substituting Eq. (B2) into Eq. (Bl), we get 

5(n ,n 2 )« ^(2 + 3 ? y + y 2 )(/o + / 2 ) + ^(4-y 2 )A, (S3) 

with y = 1/h. Therefore, the summation over the entire 
range [no,n 2 ] can be obtained approximately with only 
three values /o, /i and / 2 given. At y — > 0, Eq. (B3) 
reduces to the Simpson rule. With the approximation 
(B2), we even can carry out a summation over a part of 
the range [no,n 2 ]. For more general uses, for n < m < 
n 2 , we have 



S(no,m)^A.fo + Bf 1 + Cf 2 (S4) 



with 



A= h(y + z)[l-3z/4 + z{y + 2z)/12] 

B= h(y + z)z[l-(y + 2z)/6] 

C= -h(y + z)z[l-(y + 2z)/3]/4 

and z = (m — n )/h. 

Now, we consider the summation 5(1, oo). When f(n) 
deceases fast at n — > oo, 5(1, oo) can be obtained ap- 
proximately over a finite range with the cutoff number 
sufficiently large. We may divide this range into sev- 
eral blocks within each of which f(x) can be regarded as 
smooth function and thereby the above algorithm can be 
applied. At most cases, f(n) may vary fast at small n. 
Therefore, the stride h should be shorter at smaller n. 
Here, we introduce a point-selection program. Consider 
L successively connected blocks. The selected points di- 
vide each block into M — 1(> 2) equal-spaced segments; 
each block contains M points. The stride (the length of 
the segment) in the Z-th block is hi = h 1-1 with h a con- 
stant integer number. By such a program, the number of 
the j-th point in the Z-th block is 

,■ , M - i i M-h ,„„. 

»w] = C? - 1 + -j— j-)>* --^rr< (S5) 

for j = 1, 2, • • • , M and I = 1, 2, • • • , L. The cutoff num- 
ber is N c = ri[M,L]- By repeatedly using the above sum- 
mation scheme, one can get approximately 



5(1, oo) w ^2w p f(n p ), 



(B6) 



f(n) w /(n ) + ci(n - n ) + c 2 (n - n ) 



(B2) 



where p runs over the L(M — 1) + 1 selected points, and 
w p is the weight at point p = [j, I] ; note that because of 
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n [M,i] — such of these points should be counted 

once in the summation. The point-selection program 
uniquely determines the weights. If M > 5 is an odd 
number, applying the above summation scheme, we get 

W[i,i] = 1 

w W] = {Ah 1 - 1 - h 1 - 1 )/?,, for j = 2, 4, • • • , M — 1 
w m = {2k 1 - 1 +h 1 - l )/3, for j = 3,5,- • • , M - 2 
w [M ,i] = (2/i + 3 + fe- 1 )/6 

W[M,i]= = + /i _1 )/3 + (1 + /i)/i"76, 

for Z = 2, • • • , L - 1 
w [m ,l]= /i 1 - 1 ^ + 1/2 + /t 1 - i /6. 

To justify the above summation scheme, we here give 
an example. Consider the summation 



.?=i 



4 f 



1 



1. 



(57) 



Applying the above scheme with h =2, L = 7, and M 
9, we have 



2uv, 



4n? - 1 



= 0.99952. 



(B8) 



The relative error is (Sum - S)/S= -4.8 x 1CT 4 . Note 
that the terms under the summation in Eq. (B7) de- 
creases as l/2j 2 at large j. If one uses the known result 

00 

Er a = **/6, 

the accuracy of the summation can be improved much 
better. Instead of Eq. (B8), we calculate the following 
summation 



Sum = ^w p ( j Jl_ _ _L) 4 ^ /12 . 



(S9) 



With such an arrangement, the value of the brackets 
in the p-summation decreases as 0(n~ 4 ) at large n p . 
This summation gives a very accurate result Sum = 
1.00000015, with a small relative error 1.5 x 10~ 7 only. 

In some cases, with the given points n^/j selected in 
advance, we need to calculate the summations, 5(1, n), 
with ny i ( c ] < n < n\j +i,i c ] and i c < L. In these cases, 
because the summation over a range needs three points 
at least, the terminate number n c = n\j c ,i c ] is then deter- 
mined as follows 



3c = 



3, if jo = l 

jo + 1, otherwise. 



Since the summation 5(1, n) can be expressed as 
5(1,71) = 5(l,n[i i ; o i — 1) + S(nn ,i c ],n), the weights at 
the points n p < ntijj as tabulated above is unchanged, 
while at the points n^y, for j = l,---,j c , w p should 
be reevaluated according to the scheme as given by Eq. 
(B4) 

For testing the accuracy of the scheme for summations 
of finite terms, we consider the following example 



^ = E 



2n 



2n + l 



(B10) 



With the [h, L, M] — [2, 7, 9] program, the summation is 
approximated as 



2w n 

Sum = 



An 2 v - 1 



(£11) 



where we use the superscript n indicating the n- 
dependence of the weights. The results for Sum/s n are 
given in Table I. The relative error R.E. = (Sum — 
S n )/S n is less than 10~ 4 . 

Finally, we give the expression for the susceptibility 
X- Since it is even for Z m , we only consider the case of 
m > 0. In real space, it is given by 

00 

X (r, Z m ) = T £ G(r, z n )G(r, z n + Z m ) 

n— — 00 



= T{2 Y, G(r,z n )G(r,z n + Z m ) 

n=l 

[m/2] 

+ 2 £ G(r,z n )G(r,z n -Z m ) 

71=1 

+ G(r,z n )G(r,-z n )\^ m ™°$} 
= T{2^w p G(r,z np )G(r,z np + Zm) 

V 

+ 2 £ w [ p m/2] G(r, z np )G(r, z np - Z m ) 
v 

+ G(r,z n )G(r,-z n )\^ m ^} 
+ 5 X (r,Z m ) 

where [m/2] is the integer part of m/2, and the last term 
is given by 

00 

5 X (r,Z m )=2T G(r,z n )G(r,z n + Z m ). 

n=N c + l 

Because of G(r, z n ) — > 5 r0 /z n at n — > 00, we have 

iV c +ro 



<*X(r, ^m) = < 



E 



mir 2 T 2n-l' 
n=A r c + l 



if m > 



Mra [zd _ V - 1 if m — 



n=l 
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Note that G(r, z n — Z m ) = G*(r, Z m — z n ). Therefore, the 
Green's function at negative Matsubara frequency can be 
determined from its complex counterpart at the positive 
frequency. For a given point-selection program, G(r, z Up ) 
is known. Those values at Z m ± z n can be evaluated 
by interpolation, or G(r,z n + Z m ) w S r0 /(z n + Z m ) if 
n + m > N c . 

APPENDIX C 

In this Appendix, we discuss the problem of inverse 
Fourier transform of Pq R (c[, 0). A typical result for 
PQ ff (q, 0) is shown in Fig. 11. This function behaves 
as P cff (q, 0) oc 1/ Va 2 + Q 2 (with q 2 = q\ + q 2 y ) at q -» 0. 
The constant a vanishes at T < T c . Even at T > T c but 
close to T c , a is very small. Physically, it means that the 
pairing fluctuation is defined in a long range in real space. 
Especially, at T < T Cl the range is infinite. Therefore, its 
primary form is not suitable for numerical inverse Fourier 
transform on a finite lattice. 




FIG. 11. Function P§ s {q) in unit of 2t at Z m = 0, 
8 = 0.125, U/t = 5, and T/t = 0.0124. 

The function P o eff (q,0) can be divided into the 'singu- 
lar' part c/ ' \J a 2 + q 2 (with c a constant) and the regular 
one. There is no problem in the inverse Fourier trans- 
form for the latter one. For the 'singular' part, the task 
is to calculate the integral 



F(jx,jy) = dq x dq y 
Jo Jo 



cos(q x j x ) cos(q y j y ) 



(CI) 



where j x and j y are the coordinates of a lattice site. Re- 



peatedly integrating by part, we get 

XI jy) jy Jo dq x J Q n dqyfx (q) cos(q x j x ) cos(q y j y ) 
+ (-±y v fo dq x f2(q x ) cos(q x j x ) 
Jx So dq x f 3 (q x ) sin(q x j x ) 

-(-lpVsW, 



with 



(C2) 



h (?) = Qy Mq y + \/a 2 +~q 2 ) - ^ a 2 + q 2 



fi{q x )= In(7r + v/a 2 + it 2 + q 2 ) 

h{(lx)= 9x(ln yja 2 + q 2 x - 1) + a atan(^). 

By this way, all /'s are regular functions. 

However, because the large factors j 2 and j x at long 
distances, we need to carry out the inverse Fourier trans- 
forms in Eq. (C2) with high accuracy. The numerical 
method of the fast Fourier transforms amounts to ap- 
plying the trapezoidal rule to the integral in Eq. (C2). 
One may use a very dense mesh in the Brillouin zone for 
the transforms. But, this is uneconomical in the present 
numerical process since such transforms need to be re- 
peatedly performed. In fact, errors in the numerical in- 
tegration stem mainly from the rapid oscillation behavior 
in the integrand. Here, we present our scheme for these 
integrals in Eq. (C2). Essentially, we need to deal with 
the following integrals, 



F c {j) = / dqf{q)cos(qj) 



F s (j) 



f dqf{q) 
Jo 



sm(qj) 



(C3) 



(C4) 



where f(q) is a regular and smooth function over (0, n). 
Firstly, we consider the simple case, 



ri3 

I(qi,Q3)= / dqf(q)cos(qj) (C5) 
Jqi 



where (91,^3) is a small range. The middle point is q 2 - 
Within this range, f(q) can be expressed as 



f(q) w ai + a 2 (q - qi) + a 3 (q - q^ 2 



(C6) 



The constants a's are determined by the values fj = 
fiQj), 

ai = fi 

a 2 = (-3/i+4/ 2 -/ 3 )/2/i 
a 3 = (/! - 2/ 2 + f 3 )/2h 2 , 
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where h = q 2 — q\. Now, repeatedly integrating Eq. (C5) 
by part and use Eq. (C6), we obtain 

^"(51,93) = Ci cos(q 1 j)+C 3 cos(q 3 j)+S 1 sin(gij)+S , 3 sin(g 3 j) 

(C7) 

with 

Ci = (3/i - 4/ 2 + / 3 )fe/2x 2 
C 3 = (/i-4/ 2 + 3/ 3 )V2x 2 
Si = [/ 3 - 2/ 2 - (z 2 - l)/i]ft/x 3 

S 3 = [(z 2 - l)/ 3 + 2/ 2 - /i]ft/x 3 

and x = j ft. It is expectable that the result given by Eq. 
(C7) is more accurate than the trapezoidal rule. 

Now, dividing the range (0, it) into 2M equal spaced 
pieces, we have 

M 

F c (j) = J2 I (^k-uq2k+l) (C8) 

k=l 

with q k — (k — l)ft and (fcM+i = tt- Using the result as 
given by Eq. (C7), we get 

M 

F c (j) = 2 Wl (x)[2 £ f 2k -i cos(g 2 fc-ij) + /1 + (-IK/2M+1' 

fc=2 
M 

+4w 2 (x) J2 / 2 fcCos(g 2 fcj) 
fc=l 

(C9) 

where the functions w\ and w 2 are given by 
Mx)= 3 | T [3-^M+cos(2 a; )], 

w 2 (z) = ^(^-cosx). 

Define a new discrete function 

g k = (-l) k fk, for fc = 1,---,2M + 1. 

With this definition, equation (C9) can be rewritten as 

F c (j) = Wl {x){C 3 [f] C,[g}} + w 2 (x){Cj[f} + C 3 [g}} 

(CIO) 

where Cj [/] is the cosine Fourier transform of function / 
defined by 

JV-l 

Cjifl = 2^fk cos(q k j) + h + (-l)Viv 

k=2 

with TV = 2M + 1. Therefore, the function F c (j) can be 
evaluated by the FFT via Eq. (CIO). 
Similarly, one can get 

F s (j) = w 1 (x){S j \f\ S 3 [g}} + Mx){SAf] + s 3 [g}} 

+ |[l + ^- i ^ £l ][/i-(-l) j /iv] 

(Cll) 



where Sj[f] is the sine Fourier transform of function / 
defined by 

JV-l 

S j [f\ = 2Y,fkMQkj). 

k=2 

APPENDIX D 

In this Appendix, we rewrite the eigen-equation (6) 
in a form more convenient for the numerical calculation. 
We intend to solve the equation in real space in order 
to get rid of the prohibitive storage requirement for the 
coefficient matrix in momentum space. 

We can apply the frequency-summation scheme just 
developed in Appendix B to the present case, so reducing 
the memory size. However, to solve the eigen-equation 
in momentum space still requires tremendous memory 
size. In some cases, fortunately, the pairing function is 
short-ranged in real space. We therefore solve the equa- 
tion in real space. To transform Eq. (6) into real space, 
one needs to maintain the matrix of the coefficients to be 
symmetrical. In the follows, we present the transforma- 
tion procedure. 

(A) Define functions /(k, n p ) and ip(k, n p ) as 

/(k, n p ) = ^Tw p G(k,z np )G(-k,-z np ), (Dl) 

V>(k,n p ) = f(k,n p )(p(k,z np ), (D2) 

where w p is the weight at frequency z Up as introduced in 
Appendix B. In real space, Eq. (6) is transformed to 

y^A(r,n p ;r',n p ,)il)(r',n pl ) =A^(r,n p ) (D3) 

r'p' 

with 

A(r,n p ;r',n p >) = ^/(r - R, n p )W(R, n p , ry)/(R - r',ry) 

R 

(DA) 

and 

W(R,n p ,n p ,) = Vp(R,z„ p - z Upl ) + V P (R,z np + z„ pl ). 

(D5) 

(B) Further more, because of the lattice symmetry, we 
need only to consider the lattice sites [r] of < r y < r x . 
Those lattice sites of r x = r y are excluded since the d- 
wave pairing is under consideration. Define 

y{r,n)=2^-i>(r,n) (D6) 
F(r,R,n) = -=L=Vs 9 /( 5 r-R,n) (D7) 
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where d r — 1 + S r o, the g-summation runs over the op- 
erations of group Ci v , s g — il is the sign factor of the 
d-wave function <j>(r, z n ) under the operation g, and gr 
denotes a site coming from r operated by g. Accordingly, 
define the new matrix 

M(r 

t Tlpi Y , Tip' ) = Y, F (*, R, n p )W(R, n p , ry )F (r', R, ry ) 



[R] 

(D8) 

where again [R]-summation runs over those lattice sites 
of < R y < R x . By so doing, the eigen equation reads, 

Y M ( r > n p; r ': n P')y( r \ n P') = A J/( r , %>)■ (- 09 ) 
[p']p' 



3.00 
2.50 
2.00 
1.50 
1.00 
0.50 
0.00 
-0.50 



(It) <H.r x , Zl ) 
20(20" 2 /(r,l) 



2 



» » « 



10 



12 
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FIG. 12. Sketch of the lattice sites used for solving Eq. 
(D9). Solid circles show the reduced region. 

A sketch of the lattice sites is shown in Fig. 12. The 
reduced region [r] is taken as the shaded sites, which 
is suitable for the description of d-wave pairing. In our 
numerical calculation, the total number of sites of [r] is 
N r = 49. The normalization condition for y(r, n p ) is 



[r]p 



r,n„ 



1. 



(DW) 



By this condition, the coupling constant v is given as v = 
A/2. Shown in Fig. 13 are the typical results for f(r x , 1) 
and (p(r x ,z\). Clearly, they are short-ranged functions. 
In passing, we compare the sizes of matrices of coeffi- 
cients required for solving the eigen-equation respectively 
in real and momentum spaces. The dimension of matrix 
M in Eq. (D8) is (N r M ) x (N r M ) = 2793 x 2793, where 
Af = 57 is the number of selected Matsubara frequency. 
However, in momentum space with a 128 x 128 mesh, 
even making use of the lattice symmetry, the dimension is 
(N k M )x(N k M ) = 118560x118560, where N k = 32x65 
is the number of momentum k with < k y < k x < n. 



FIG. 13. Functions f(r x ,l) and 4>(r x ,zi) at r y = 0, 
5 = 0.125, U/t = 5, and T/t = 0.0124. The dashed lines 
are for the eyes. 
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100. 
500. 



1.000013 
1.000013 



0.000013 
0.000013 



TABLE I. Ratio of Sum given by Eq. (All) and 
S n = 2n/(2n + 1) at various n. R.E. represents the relative 



error. 



Sum/S n R.E. 



5. 1.000000 

10. 1.000047 0.000047 

50. 1.000012 0.000012 



17 



